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1. 


ABSTRACT 


Neck  pain  is  pervasive  problems  in  military  population,  especially  in  those  working  in 
vibrating  environments.  Previous  studies  show  neck  pain  is  strongly  associated  with  degeneration 
of  Intervertebral  Disc  (IVD),  which  is  commonly  caused  by  repetitive  loading  and  aging.  To 
reduce  the  risk  of  cervical  spine  disease,  there  is  a  need  to  measure  the  effect  by  helmet,  equipment 
and  seating.  However,  in-vivo  displacement  and  loading  condition  of  cervical  spine  are  difficult  to 
measure  during  operation.  Clinical  ultrasound  (US)  is  being  explored  as  a  tool  to  image  motion  of 
the  cervical  spine,  specifically  vertebral  kinematics  and  intervertebral  disc  deformation.  A  dual 
US  imaging  system  is  being  developed  to  measure  3D  motion  of  contiguous  cervical  vertebrae  and 
IVD  strain  (Aheight  of  IVD/original  IVD  height)  of  intervening  functional  spine  units  (FSU).  This 
system  was  validated  ex-vivo  using  cadaveric  C-spines  mounted  in  a  servo-hydraulic  material 
testing  machine  by  comparing  dynamic  US  measurement  to  direct  measurements  using  the  linear 
variable  differential  transfonner  (LVDT)  of  the  hydraulic  test  system.  Over  the  past  year  of  support, 
2014-2015,  we  improved  the  perfonnance  of  hardware  and  software  for  the  dual  ultrasound  system. 
The  capability  of  the  dual  US  to  measure  C-spine  elastic  and  viscoelastic  properties  in-vivo  under 
conditions  simulating  vibrations  over  a  range  of  frequencies  and  amplitudes  is  currently  being 
tested.  Finite  element  (FE)  models  of  the  human  cervical  spine  have  been  developed  to  examine 
the  dynamic  response  of  the  cervical  spine  functional  spine  units  in  response  to  forces  and 
moments  applied  to  the  head  and  neck.. 
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2.  INTRODUCTION 


In  the  year  2013,  we  developed  a  unique  dual  US  system  that  can  non-invasively  measure 
IVD  deformation  and  mechanical  compliance  ex-vivo,  and  provide  real-time  images  of  IVDs  and 
dynamic  vertebral  motion  in-vivo  during  simulated  tasks.  A  transfer  function  that  modeled  the  non¬ 
linear  displacement  of  human  cadaveric  cervical  spine  functional  spine  units  was  developed  to 
account  for  the  differential  elasticity  during  compression  or  tension. 

The  project  focus  of  2014  was  to  apply  the  dual  US  in-vivo  and  model  analytically  the 
elastic  and  viscoelastic  behavior  of  the  cervical  spine  in  response  to  cyclic  loading.  Compared  to 
using  the  US  system  to  image  the  C-spine  ex-vivo ,  using  the  dual  US  system  to  image  the  C-spine 
in-vivo  was  much  more  challenging,  affected  by  inadvertent  human  motion  and  inconsistent 
muscle  activity  that  served  to  dampen  the  system.  2014  accomplishments: 

1.  Improved  software  for  dual  US  system  by  developing  a  mathematical  algorithm  that  corrects 
signal  drift  during  tracking  of  vertebrae  during  long  duration  ultrasound  imaging.  This 
significantly  increased  tracking  accuracy  for  in-vivo  US  imaging; 

2.  Developed  system  to  apply  in-vivo  static  and  dynamic  cyclic  loads  to  the  cervical  spine  at 
known  forces,  displacements  and  frequencies  allowing  derivation  of  the  transfer  function  for 
C-spine  FSU’s  C4-5; 

3.  Established  a  robust  protocol  for  ex-vivo  FSU/IVD  fatigue  testing.  Elastic,  viscoelastic  and 
failure  properties  of  FSU  determined  for  C-spine  FSUs  that  were  incorporated  into  material 
properties  of  C-spine  FEM  used  to  simulate  cyclic  loading  and  fatigue  behavior  of  C-spine. 
The  biomechanics  properties  (compliance  or  damping  coefficient)  of  the  IVD  have  been 
correlated  to  clinical  measures  of  FSU  health  and  function  measured  by  MRI  (Pfirrmann  grade) 
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3.  PROGRESS 


3.1  Ultrasound  System  Development  and  Validation  (Continuation  of  Task  1  and  2  in 
Aim  1) 

A  blockmatching  algorithm  was  used  to  compare  two  consecutive  US  images  and  map 
the  lateral  (i.e.  perpendicular  to  U/S  beam)  and  axial  (parallel  to  U/S  beam)  displacement  of 
pixels  fonning  the  image  profiles  of  the  vertebrae.  This  algorithm  was  previously  validated  using 
ex-vivo  cadaver  cervical  spines.  However  during  real  time  in-vivo  tracking  of  cervical  vertebrae, 
new  challenges  emerged  that  did  not  occur  in  the  cadaver  test:  it  was  difficult  to  maintain  the 
dual  ultrasound  probes  in  the  standardized  fixed  position  in  direct  acoustic  contiguity  with  the 
vertebrae  during  dynamic  flexion/ex  tension  motion.  The  previous  methods  failed  to  correct 
signal  drift,  which  was  accumulated  during  long  duration  US  tracking.  Accuracy  of  FSU 
kinematics  especially  during  high  frequency  motion  was  limited. 

Based  on  the  work  of  Rahimi  et  al.  2001,  we  developed  a  mathematical  model  suitable  to 
correct  signal  drift  during  tracking  of  vertebrae  during  long  duration  ultrasound  imaging.  This 
method  relaxes  the  restriction  that  only  temporally  adjacent  ultrasound  frames  are  compared 
(known  as  frame-to-frame  tracking)  to  calculate  the  relative  interval  motion  of  the  vertebra.  It 
allows  high-quality,  non- sequential  measurements  to  be  used  in-lieu  of  poor  quality,  sequential 
measurements.  When  estimating  the  position  of  vertebrae,  the  motion  tracking  method  takes  into 
account  the  change  in  position  between  the  current  frame  and  the  previous  frame,  as  well  as  other 
previous  spatially  related  image  frames  tracking  the  position  of  the  vertebra.  This  expanded  set  of 
high  quality  US  images  is  now  used  to  analyze  the  ultrasound  radio  frequency  (RF)  data. 

The  accuracy  and  precision  of  US  measurements  of  vertebral  motion  were  validated  using 
a  PMMA  phantom  to  track  rigid  body  motion  over  a  range  of  frequencies.  A  PMMA  block  was 
secured  to  the  actuator  piston  of  a  servohydraulic  materials  testing  machine  to  which  1mm 
sinusoidal  displacements  were  applied  over  a  range  of  frequencies  for  20  seconds  (1000  frames). 
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The  block  was  fully  immersed  in  a  water  bath  to  allow  acoustic  coupling  to  the  US  probe.  The 
applied  displacement  was  compared  between  direct  measurements  from  the  linear  voltage 
differential  transformer  (LVDT)  in  line  with  the  actuator  and  that  derived  using  a  single  US  device 
and  the  motion  tracking  software  (Figure  1). 


Figure  1:  Ultrasound  imaging  of  PM M A  phantom  secured  to  actuator  piston  of  materials  testing  machine. 

Measurements  were  compared  between  derived  displacements  using  tracking  algorithm  and  LVDT. 

Six  different  Regions  of  Interest  (ROIs)  on  the  phantom  were  chosen  to  evaluate  the 
precision  of  the  tracking  algorithm.  For  applied  frequencies  1-8  Hz,  the  mean  absolute  error  (MAE) 
between  translational  motion  measured  using  the  US  tracking  system  and  the  LVDT  ranged  from 
0.0227  to  0.0508  mm  (Figure  2).  The  overall  average  MAE  of  US  measurements  were  0.041  mm. 
No  drift  was  observed  during  the  course  of  scan  (1000  frames).  Accuracy  was  not  significantly 
affected  by  motion  frequency.  Assuming  typical  IVD  deformation  of  0.8mm  under  1 50N  load,  this 
corresponds  to  a  6.3%  precision  error  in  the  US  measured  strain. 

Software  using  the  new  tracking  algorithm  has  been  integrated  into  dual  US  package.  The 
software  has  graphic  user  interface  (Figure  3),  which  allows  the  operator  to  analyze  the  relative 
displacements  of  cervical  vertebrae  and  the  deformation  of  intervertebral  discs  (IVD).  The 
program  is  written  on  C++  under  open  source  framework  and  the  executive  files  are  compatible 
with  different  systems  (tested  in  PC  and  Mac),  which  makes  the  dual  ultrasound  system  (US) 
usable  as  a  stand-alone  system  to  acquire  cervical  spine  motion  data  to  calculate  kinematics  and 
IVD  deformation. 
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Figure  2:  The  Mean  Absolute  Error  for  US  vs.  LVDT  measurements  of  phantom  motion  plotted  as  function 
of  applied  frequency. 
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Figure  3:  Graphic  User  Interface(GUI)  of  tracking  software  which  allows  users  to  analyze  displacement  of 


objects.  Upper  left  plot  is  the  displacement  curve.  Upper  middle  window  plots  the  displacement  map 
between  frame  pairs.  Since  we  are  applying  this  program  to  track  a  rigid  phantom,  all  points  inside 
phantom  show  the  same  displacement.  B-mode  image  of  phantom  is  shown  lower  middle,  with 
measurements/image  details  shown  in  lower  right. 
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3.2  Ex-vivo  Hydration  vs  IVD  Compliance  and  Cervical  Spine  Health  and  Integrity 
(Task  4  and  6) 

Methodology  to  measure  the  mechanical  properties  of  single  and  contiguous  FSUs  using 
intact  cadaver  specimens  C2-C7  was  developed.  Previous  ex-vivo  experiments  were  conducted  on 
isolated  cervical  spine  FSU’s  which  fail  to  replicate  physiological  conditions,  especially  the 
lordotic  curvature  of  C-spine. 

Contiguous  segments  C2-C7  from  an  84-year-old  male  were  used  to  develop  the  protocol. 
The  posterior  skin  was  retained  to  provide  acoustic  coupling  for  the  dual  US  probes.  The  trachea 
and  anterior  skin  were  removed,  but  paracervical  muscles  and  ligament  were  retained  to  maintain 
the  tensile  and  compressive  properties  of  the  connective  soft  tissue  envelope.  The  C2  and  C7 
vertebral  bodies  were  potted  in  liquid  plastic.  Using  minimally  invasive  technique,  uniaxial,  5  mm 
diam.  titanium  pedicle  screws  were  placed  into  the  vertebral  bodies,  just  anterior  to  the  transverse 
process,  under  fluoroscopic  guidance.  To  maintain  the  cervical  lordosis,  a  slightly  bent,  5  mm 
diameter  stainless  steel  interconnecting  rod  was  secured  in  the  tulip  head  connectors  of  contiguous 
pairs  of  pedicles  screws,  rigidly  secured  with  set  screws.  Our  method  was  slightly  modified  from 
pedicle  screw  placement  protocol.  The  pedicle  screws  were  placed  in  the  vertebral  bodies  in  an 
Anterior-Posterior  (A-P)  position  to  better  control  (limit)  flexi on/extension  motion  of 
“instrumented”  FSU’s  during  cyclic  mechanical  testing.  By  using  variable  lengths  of 
interconnecting  rods  secured  to  different  pedicle  screw  levels,  different  combinations  of  FSU’s 
could  be  fixed  or  mobilized.  The  specimen  was  mounted  in  material  testing  system  (Instron  8511, 
Norwood,  MA)  and  loaded  from  20  lbs.  compression  to  20  lbs.  tension.  Fluoroscopic  images 
captured  the  position  of  the  vertebral  bodies  during  different  loading  conditions  (Figure  4).  Dual 
ultrasound  transducers  were  positioned  anteriorly  and  posteriorly  to  image  the  IVD  and  vertebral 
bodies  anteriorly  and  spinous  processes  posteriorly.  There  was  no  movement  between  the  fixed 


9 


FSUs  during  compression  and  tension  (Figure  5),  while  the  “empty”  pedicle  screws  did  not 


interfere  with  the  motion  of  the  mobile  FSUs. 


Figure  4:  Cadaver  testing  with  dual  US  and  C-arm  X-ray. 


Figure  5:  X-ray  images  of  20lbs.  of  tension  (left)  and  20lbs.  compression  (right).  Black  bars  are  US 
transducers.  Relative  motion  occurred  between  mobile  FSU  levels  (red  bar  length  change  =  0.1435mm), 
while  no  motion  occurred  between  fixed  FSU  levels  (yellow  bar  length  change  =  0. 029mm,  image  resolution 
of  X-ray  is  0.3mm). 


The  stiffness  and  damping  coefficients  of  the  C4-C5  FSU  (most  common  level  of 


degenerative  disease)  were  correlated  with  qualitative  grading  (Pfirrmann  scale)  of  the  IVD  on 


MRI  images  (Figure  6).  The  ability  of  the  dual  US  system  to  measure  the  biomechanical  properties 


of  the  FSU  was  further  validated  by  tracking  the  deformation  of  the  C4-C5  FSU  in  real  time  during 


a  standard  creep  test.  Five  fresh  frozen,  cadaveric  human  cervical  spines,  C2-T1,  ages  54  -  67  years 


(Medcure,  Portland,  OR)  were  defrosted,  the  trachea,  esophagus,  and  skin  were  removed.  The 


surrounding  para-cervical  muscles  and  adipose  tissue  were  retained.  Each  spine  was  subjected  to 
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150N  compressive  creep  test.  Using  minimally  invasive  technique,  uniaxial,  5  mm  diam.  titanium 
pedicle  screws  were  placed  into  the  vertebral  bodies  under  fluoroscopic  guidance.  Steel  rods  were 
used  to  fix  C2-C4  FSUs  and  C5-T1  FSUs.  The  deformation  of  the  C4-C5  IVD  was  measured  by 
the  dual  US  system  and  correlated  with  the  real-time  applied  load. 
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Figure  6:  B-mode  images  of  Ultrasound  are  compared  to  MR!  sagittal  image  of  cadaveric  C-spine. 
Ultraound  provides  a  magnified  window  of  Functional  Spinal  Units  in  cadaver.  One  US  acquired  anterior 
vertebrae  body  and  the  other  imaged  posterior  spinal  processes. 

Stiffness  and  damping  coefficients  of  C4-C5  FSU  were  derived  assuming  a  standard  Voigt 
model  comprised  of  an  elastic  spring  and  a  viscous  damper  connected  in  parallel.  Multivariate 
Linear  Regression  was  used  to  calculate  the  stiffness  of  the  elastic  spring  and  the  damping 
coefficient  of  viscous  damper.  These  derived  mechanical  properties  were  compared  to  the 
structural  integrity  of  the  IVD  measured  on  T2  weighted  MRI  images  of  each  FSU  based  on  the 
Pfirrmann  grading  system,  which  is  the  clinical  standard  for  evaluating  the  “health”  of  the  IVD. 
Grade  1  is  healthy  hydrated  disc  that  shows  homogeneous  brightness  on  T2  MRI  image  while 
black,  collapsed,  dehydrated  IVD  on  T2  MRI  image  is  graded  5. 

The  age  of  the  cadaveric  C-spine  specimens’  ranged  from  50-70  years,  and  the  Pfirrmann 
grades  ranged  from  2  to  4.  For  compressive  creep,  the  stiffness  of  the  C4-5FSU  correlated  better 
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with  the  Pfirrmann  Grade  than  the  damping  coefficient  (Table  1).  Younger  specimens  with  better 
IVD  hydration  and  structural  integrity  tended  to  be  more  compliant  compared  to  older  specimens 
with  higher  Pfirmann  grades. 


Table  1  Material  properties  of  C4-C5  FSU  measured  by  dual  ultrasound  system 


Specimen# 

1 

2 

3 

4 

5 

Stiffness  (N/mm) 

117.3 

203.3 

258.3 

264.2 

312.5 

Damp  Coefficient  (N*s/mm) 

1420 

1497 

1045 

1534 

1110 

Pfirrmann  Grade  of  IVD 

2 

3 

3 

3 

4 

Age  of  Specimen  (years) 

54 

58 

63 

58 

67 

3.3  In-vivo  Ultrasound  Test  in  Simulated  Environment  (Task  7) 

In  order  to  derive  a  transfer  function  for  the  mechanical  behavior  of  the  cervical  spine  in- 
vivo  in  a  consistent  and  reproducible  way,  we  are  developing  a  system  that  applies  cyclic 
distraction/compression  loads  to  the  neck  at  known  amplitudes  and  frequencies  representative  of 
those  mounted  troops  are  exposed  to  during  military  operations.  Two  different  test  systems  have 
been  explored,  one  that  actively  applies  cyclic  distraction  to  the  head/neck,  similar  to  cervical 
traction  used  to  treat  cervical  injuries  clinically  and  another  system  that  applies  a  cyclic 
compressive  load  to  the  head/neck  via  an  actuator  lifting  the  seat  in  which  the  subject  is  sitting, 
more  representative  of  the  scenario  for  mounted  troops.  The  distraction  based  system  consists  of 
a  motorcycle  helmet  worn  by  the  subject  with  an  attached  weight  (101b.)  that  applies  a  static 
compressive  load  to  the  cervical  spine.  The  weighted  helmet  is  connected  by  a  cable  to  a  solenoid 
powered  actuator  that  “unloads”  the  spine  by  applying  a  variable  traction  load  (0  to  20  lbs)  to  the 
head  and  neck  in  a  cyclic  fashion  over  a  range  of  frequencies  (Figure  7,  Figure  8).  A  load  cell 
interposed  between  the  weight  and  the  helmet  measures  the  resultant  force  applied  to  the  neck  in 
real  time  (iLoad  Mini,  Loadstarsensor,  Fremont,  CA). 
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Figure  7:  Schematic  of  cyclic  load  application  system  to  head/neck  of  subject  while  sitting  or  standing  in 
load  frame. 


Figure  8  Solenoid  powered  traction  system  for  applying  cyclic  loads  in-vivo 


To  ensure  that  the  loading  axis  of  the  head  and  neck  were  collinear  with  the  applied  uniaxial 
force,  a  5  point  latch  &  link  restraint  system  was  used  to  fix  the  torso  of  the  sitting  subject  to  the 
load  frame  along  the  line  of  action  of  the  cable.  We  tested  the  solenoid  based  traction  system  on 
ourselves  by  perfonning  static  and  dynamic  analysis  of  the  compliance  of  C4-C5. 
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Figure  9:  Simplified  rigid  body  representation  of  cervical  spine  and  coordinate  system.  The  anterior  profile 
of  the  vertebral  body  and  IVD  is  imaged  by  the  anteriorly  positioned  US  probe,  while  the  facet  joint  (cross- 
hatched  region)  and  spinous  process  is  imaged  by  the  posteriorly  positioned  US  probe  (~60  ° relative  to  the 
anteriorly  positioned  probe) 

The  geometry  of  the  vertebrae  is  assumed  to  remain  constant  during  the  course  of  the  test. 
The  surface  contours  of  the  cervical  vertebrae  imaged  by  the  dual  ultrasound  probes  can  be 
depicted  as  a  simplified  2D  rigid-body  (Figure  9).  A  global  Cartesian  coordinate  system  is  fixed 
relative  to  the  torso.  The  anterior  and  posterior  motion  of  the  C4-C5  FSU  is  represented  by  two 
vectors  (vertebral  body  endplate  motion)  and  i>  (posterior  element  motion). The  deformation  of 
the  C4-C5  IVD  is  estimated  from  the  average  displacement  of  these  two  vectors,  while  the 
flexion/extension  angle  of  C4-C5  is  estimated  by  calculating  their  vector  cross-product: 


Displacment  —  ^TA^rp^ 


.  ,  ^Ay  ^yO 

Angle  —  arctan - 

T~Ax  Txq 


+  arctan 


rAx>  rAy>  rPx>  r py  are  the  components  of  r^r^in  the  superior-inferior  (X)  and  anterior- 
posterior  (Y)  directions,  respectively  in  the  loaded  configuration  and  rx0  and  ry0  are  the 
components  of  these  vectors  in  the  unloaded  configuration. 

Static  Testing:  A  static  force  from  -10  lbs.  (compression)  to  +15  lbs.  (traction)  was  applied 
incrementally  to  the  cervical  spine  using  our  controlled  traction  system  (Figure  10).  The  dual  US 
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system  measured  the  resultant  C4-C5  FSU’s  displacement  and  angular  motion  (Figure  10).  The 
static  compliance  was  derived  from  the  slope  of  the  displacement  vs.  load  curve  deduced  by  linear 
regression  (Figure  11). 


Figure  10:  Sequential  step  loads  applied  to  the  head  and  neck  witht  resultant  C4-C5  displacement  and 
flexion/extension  angle  measured  by  US.  The  morphology  of  load-time  curve  and  displacement-time  curve 
demonstrated  that  the  IVD  displacement  at  C4-C5  correlated  with  the  applied  axial  force.  The  resultant 
angular  deformation  implies  that  the  applied  axial  force  also  exerts  a  moment  on  the  cervical  spine,  a 
consequence  of  either  the  natural  cervical  spine  lordosis,  misalignment  between  the  line  of  action  of  the 
applied  load  and  the  neutral  axis  of  the  cervical  spine  or  voluntary  muscle  contraction  by  the  subject. 


Figure  11:  Elastic  compliance  of  C4-C5  FSU  using  solenoid  based  traction  device  to  apply  static  incremental 
loads  R2=  0.97(red,  25y/o  subject),  0.76  (blue,  56y/o  subject).  FSU  tended  to  be  more  compliant  for  the 
younger  human  subject  (0.00589mm/N)  compared  to  older  subject  (0.00265mm/N). 

Dynamic  Testing:  The  feasibility  of  using  US  to  measure  the  dynamic  mechanical  properties  of 
cervical  spine  FSUs  was  explored  by  cyclically  loading  the  neck  between  0  to  +20  lbs.  The  applied 
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load  was  a  square-wave  function  (Figure  12a),  the  resultant  displacement  was  highly  non-linear 
(Figure  12b)  exhibiting  two  distinct  phases  during  loading  and  unloading:  1)  there  was  little  FSU 
displacement  during  the  application  and  release  of  the  traction  load  impulse;  2)  there  was  creep 
deformation  during  the  phase  that  the  applied  load  was  held  constant.  Analysis  in  the  frequency 
domain  facilitated  derivation  of  the  dynamic  compliance,  while  the  damping  effect  was  calculated 
from  the  area  under  the  displacement  -  load  curve  which  physically  represents  the  work  lost  in 
each  cycle. 


Figure  12:  (a)  Cyclic  tensile  step  load  applied  by  solenoid  actuator;  (b)  resultant  C4-C5  FSU  deformation 
measured  by  dual  US  system. 

Dynamic  frequency  analysis  demonstrated  the  relative  power  spectra  of  applied  loads  and 
displacements  (Figure  13).  The  highest  peak  in  the  frequency  spectrum  corresponded  to  the  load 
and  displacement  at  the  applied  frequency  of  the  square  wave  and  was  used  to  fdter  out  higher 
frequency  “noise”.  Dynamic  compliance  was  derived  from  the  ratio  of  the  C4-C5  FSU  defonnation 
in  response  to  the  applied  load  at  the  applied  frequency.  The  25  y/o  subject  had  a  more  compliant 
C4-C5  FSU  than  the  56  y/o  subject  at  all  frequencies,  but  the  compliance  decreased  at  higher 
frequencies.  The  compliance  of  the  C4-C5  FSU  for  the  56  y/o  subject  changed  little  as  a  function 
of  the  applied  frequency  (Table  2).  The  damping  coefficient  was  more  affected  by  frequency  than 
age  (Table  3). 
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Figure  13:(a )  Frequency  analysis  of  applied  load  waveform  (b)  C4-C5  displacement  measured  by  dual  US 
in  frequency  domain 


Table  2  Dynamic  compliance  as  function  of  applied  frequency  and  age 


Frequency  applied 

load 

Dynamic  Compliance  mm/N 

1  Hz 

2Hz 

3Hz 

Younger  Subject 

0.0306 

0.0242 

0.0098 

Older  Subject 

0.0111 

0.0147 

0.0123 

Table  3  Work  Loss  as  Function  of  Applied  Frequency  and  Age 


Frequency  of  load 
applied  to  neck 

Area  Estimate  (Work  loss)  N  x  mm 

1  Hz 

2Hz 

3Hz 

Younger  Subject 

39.28 

64.81 

75.03 

Older  Subject 

42.95 

64.23 

127.13 

While  this  data  is  preliminary  and  limited  to  the  PI  and  his  graduate  student,  it  demonstrates 
the  promise  of  the  proposed  dynamic  actuator  to  apply  static  and  dynamic  loads  to  the  C -spine  in 
a  controlled,  parametric  fashion  while  simultaneously  measuring  the  resulting  deformation  of  C- 
spine  FSU’s  using  the  dual  US  system  thereby  allowing  calculation  of  the  mechanical  properties 
of  the  cervical  spine  in-vivo  that  dictate  its  functional  performance.  Differences  in  the  transport  of 
water  through  the  poroelastic  matrix  of  the  IVD  (which  will  be  investigated  further  using  diffusion 
based  MRI)  may  account  for  the  apparent  affect  of  age  and  applied  frequency  on  the  compliance 
of  the  FSU.  Additionally,  reflexive  contraction  of  the  paracervical  muscles  to  dampen  forces 
applied  to  the  cervical  spine  at  higher  loading  frequencies  could  also  account  for  the  observed 
decrease  in  compliance  and  increase  in  energy  dissipation  at  higher  applied  loading  frequencies. 
EMG  measurement  of  the  activation  of  the  paracervical  muscles  will  elucidate  this  further. 
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These  preliminary  data  identified  a  significant  shortcoming  of  the  solenoid  actuator  to 
apply  cyclic  distraction/compression  loads  at  frequencies  and  amplitudes  that  will  facilitate 
derivation  of  a  universal  transfer  function  for  the  mechanical  behavior  of  the  cervical  spine  FSUs 
in-vivo  in  a  consistent  and  reproducible  fashion  or  allow  the  application  of  different  load  profiles 
(sinusoidal,  square  wave,  triangle,  combinations  thereof)  that  simulate  the  forces  and  moments 
applied  to  the  head  and  neck  of  mounted  troops  during  military  maneuvers.  Therefore  we  have 
developed  a  new  system  that  uses  a  pneumatic  powered  load  actuator  that  will  allow  more 
controlled  cyclic  load  applications  (Figure  14). 

In  this  new  design,  the  force  is  applied  vertically  through  the  seat  that  the  subject  is  seated, 
transferred  up  the  spine  to  the  head  and  neck,  which  is  secured  in  a  helmet  rigidly  fixed  to  an 
adjustable  cross-bar.  At  the  point  of  attachment  of  the  helmet  to  the  cross-bar,  there  is  a  safety  clutch 
that  “breaks”  for  any  applied  load  exceeding  25  lbs  tension  or  compression  (Figure  14B).  Four  co¬ 
axial  linear  bearing  attached  to  the  seat  ensure  precise  up/down  vertical  movement  without  binding 
(Figure  14C).  The  maximum  extension  of  the  actuator  is  limited  by  a  separate  position  sensing  system. 
To  safeguard  against  the  risk  of  irregular  air  flow  to  the  actuator,  two  redundant  valves  were  added 
with  an  electronic  PSI  regulator.  Additional  safety  features  include  an  emergency  “stop”  switch  that 
can  be  activated  by  either  the  subject  or  technician.  Validation  of  the  system  and  a  full  safety 
evaluation  is  being  conducted  using  a  Crash  dummy  (Humanetics,  Plymouth,  MI),  an  established  test 
device  in  accordance  with  the  USA  Code  of  Federal  Regulations.  The  dummy  (Figure  15)  has  two  6- 
axis  load  cells  on  the  top  and  bottom  of  its  simulated  C-spine,  which  provide  real  time  load  data 
applied  to  the  neck.  Preliminary  studies  on  the  break-away  device  were  performed.  The  break-away 
device  consists  of  4  sets  of  commercially  available  wire  break-away  connecters  (DCD  Design,  British 
Columbia,  Canada),  that  break  when  the  applied  load  exceeds  the  strength  of  internal  wire.  Repetitive 
loading  studies  that  applied  hyper-extension  and  hyper-compression  to  the  dummy’s  spine 
demonstrated  that  the  device  released  at  -245. 0±9. 1  N  (~  50  pounds)  compressive  load  for  a  maximum 
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applied  flexion  moment  of  9. 1+3.4  Nm  and  291±28  N  (~  60  pounds)  for  a  maximum  applied  extension 
moment  of  14.1+8.5  Nm  (Table  4).  This  exceeds  our  design  criteria  of  25-301bs  maximum  load  to  be 
applied  to  the  spine  elements.  Therefore  customized  wires  from  manufacturer  will  be  used  to  attain 
this  safety  standard. 


Figure  14:  (A)  front  view  of  the  updated  design;  (B)  safety  device  release  when  applied  force  exceeds 
limit;  (C)  Lateral  view  of  linear  bearings  and  actuator  placement 


Figure  15:  lateral  view  of  the  dummy  in  spine  cyclic  loading  system  and  its  Cartesian  coordinates  system 
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Table  4  Force  and  Moment  Applied  to  Neck  When  Safety  Device  Breaks  away 


Extension 

Compression 

Title 

Average 

STDEV 

Avgrage 

STDEV 

Upper  Neck  (Cl) 

x-Axis  Force(N) 

66.6 

8.7 

-67.9 

6.1 

y-Axis  Force(N) 

11.2 

8.8 

-7.5 

1.3 

z-Axis  Force(N) 

272.3 

28.8 

-232.2 

13.4 

x-Axis  Moment(Nm) 

1.0 

0.8 

-0.3 

0.1 

y-Axis  Moment(Nm) 

-1.7 

3.1 

-3.3 

1.1 

z-Axis  Moment(Nm) 

0.7 

0.4 

-0.7 

0.2 

Lower  Neck  (C7) 

x-Axis  Force(N) 

61.2 

9.4 

-53.0 

9.7 

y-Axis  Force(N) 

15.0 

5.2 

-6.6 

1.0 

z-Axis  Force(N) 

290.9 

27.9 

-245.0 

9.1 

x-Axis  Moment(Nm) 

1.8 

1.1 

-0.7 

0.2 

y-Axis  Moment(Nm) 

14.1 

8.5 

-9.1 

3.4 

z-Axis  Moment(Nm) 

1.0 

1.0 

-0.5 

0.0 

Cyclic  loading  of  the  dummy  was  attempted  at  different  frequencies  (lHz-5Hz).  At 
frequencies  higher  than  2Hz  the  pneumatic  actuator  was  unable  to  apply  a  symmetric  load- 
displacement  curve;  the  rate  of  the  applied  load  to  lift  the  seat  and  dummy  was  slower  than  the  rate 
of  unloading.  This  is  a  result  of  the  difference  in  the  rate  of  pressurizing  and  depressurizing  the 
pneumatic  actuator.  To  compensate,  we  will  use  dead  weights  equivalent  to  the  weight  of  the  subject 
and  the  seat  to  “tare”  the  system,  so  that  the  pneumatic  actuator  only  applies  ±  251bs.  Cadaveric  torsos 
with  attached  cervical  spines  C2-C7  will  also  be  tested  in  this  device  prior  to  any  in-vivo  human  tests 
to  further  establish  the  capabilities  of  the  system  to  measure  cervical  spine  FSU  kinematics  and  IVD 
deformation.  These  data  will  be  compared  to  similar  tests  conducted  on  isolated  cervical  spines  using 
the  Instron  servohydraulic  test  system. 
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3.4  FE  Model  Validaition  and  Modification  (Task  11)  and  Simulate  Operation  Loading 
Conditions  (Task  12) 

Finite  element  (FE)  modeling  for  prediction  of  fatigue-related  material  property  changes 
in  cervical  spine  soft  tissues  is  ongoing.  N-code  (HBM,  Southfield,  MI),  was  used  to  detennine 
fatigue  in  the  disc  and  it  was  able  to  reproduce  fatigue  in  terms  of  number  of  cycles  for  a 
prescribed  compressive  load.  However  the  formulation  for  calculation  of  Stress  -Number  of 
Cycles  (SN)  curve  used  is  N-code  is  suitable  for  metals  than  for  soft  tissues  as  it  is  directly 
calculated  from  initial  stiffness  of  the  material,  hence  it  is  difficult  to  validate  the  results. 
Considering  the  observations,  conclusions  and  challenges  from  the  previous  exercise  a  unique 
method  of  FE  modeling  was  developed  during  the  last  quarter.  The  new  method  consists  of 
coupling  of  Matlab  (The  Mathworks  Inc.,  Natick,  MA)  and  LS  Dyna  (Livermore  Software 
Technology,  Livennore,  CA)  to  carry  out  fatigue  simulations.  This  method  was  based,  in  part,  on 
Qasim  et  al.  2012  and  Qasim  et  al.  2014,  who  adopted  a  similar  technique  using  Fortran  in 
ADINA  (ADINA  R&D  Inc.,  Watertown,  Massachusetts). 

Methodology:  The  computational  algorithm  which  was  used  to  model  material  property  changes 
in  cervical  spine  soft  tissues  related  to  repetitive  axial  compressive  loading  was  outlined  in 
Figure  16.  Fatigue  experiments  with  varying  force  amplitude  and  frequency  were  performed  on 
bone-disc-bone  segments  at  an  interval  of  10,000  cycles,  up  to  50,000  cycles.  Quantified  changes 
in  intervertebral  disc  properties  included:  cycle-by-cycle  compressive  stiffness,  intervertebral 
disc  height,  viscoelasticity,  and  tensile/compressive  stiffness.  The  experimental  results  were  used 
as  input  to  the  FE  model  as  change  in  material  properties  due  to  fatigue. 
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Figure  16:  Schematic  representation  of  vertebral  disc  segment  under  repetitive  loading  in  a  saline  bath. 

A  FE  model  of  a  single  cervical  spine  segment  was  developed.  The  IVD  model  consists  of 
annulus  fibrosus,  nucleus  pulposus,  and  superior  and  inferior  endplates.  The  IVD  model  was 
discretized  and  programed  with  element- specific  material  properties  with  a  script  in  Matlab.  This 
process  changing  of  material  property  values  to  model  fatigue-related  changes  based  on  element- 
specific  loading  conditions  after  each  cyclic  loading  procedure.  The  procedure  consists  of  first 
performing  a  computational  set  of  fatigue  cycles  in  LS  Dyna,  then  using  a  Matlab  script  to  alter 
element- specific  material  properties  of  the  soft  tissues  using  Von  Mises  stress  as  the  decision 
parameter,  followed  by  the  next  fatigue  cycle  simulation  in  LS  Dyna.  Changes  in  IVD  material 
properties  are  based  on  experimental  results  and  will  include  scaling  factors  for  the  rate  and 
magnitude  of  applied  repetitive  compressive  loads.  The  degraded  material  properties  will  be 
incorporated  for  elements  with  Von  Mises  stress  higher  than  the  average.  This  report  highlights 
experimental  results  from  cervical  spine  segment  testing  to  quantify  fatigue-related  material 
property  changes  and  discretization  of  the  finite  element  model,  along  with  the  experiments  and 
FE  fatigue  modeling  done  before. 

PART  1:  Experiments 

Experiments  performed  during  the  three  quarters  of  last  year  are  described  below.  Based 
on  the  experience  of  results  from  the  previous  set  of  experiments,  the  experimental  protocol  was 
updated.  The  experiments  can  be  considered  to  be  performed  in  three  parts  during  the  year  and  the 
protocol  adopted  during  quarter  4,  yielded  best  results  and  will  be  used  for  the  future  experiments. 
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Part  1:  The  experiments  were  performed  to  understand  the  effect  of  pure  moment  tests, 
viscoelastic  tests  in  between  fatigue  tests  and  the  inclusion  of  saline  bath.  Pure  moment  testing 
was  first  used  to  measure  the  bending  response  of  cervical  segments.  Testing  consisted  of  the 
static  application  of  pure  moments  at  0.33Nm,  0.5Nm,  INm,  1.5Nm  and  2Nm  under  flexion- 
extension  and  lateral  bending  (Wheeldon  et  al.  2006).  Following  pure  moment  testing,  the  segment 
was  mounted  on  an  anvil  of  MTS  Systems,  Eden  Prairie,  MN,  which  was  immersed  in  a  saline 
bath  maintained  at  34oC.  The  segment  was  allowed  to  acclimate  in  the  water  bath  for  1  hour. 
Stress  relaxation  testing  was  then  performed  to  quantify  viscoelasticity  of  the  segment.  The 
specimen  was  compressed  by  5%  of  the  disc  height  and  allowed  to  relax  under  the  constant 
deformation.  Stress  relaxation  testing  was  also  performed  prior  to  and  following  all  fatigue  sets. 
A  Maxwell  material  model  under  fixed  compressive  strain  was  used  to  quantify  viscoelasticity 
during  the  stress  relaxation  testing.  Segments  were  then  exposed  to  cyclic  loading  between  0  and 
-150N  for  1000,  10,000  and  20,000  cycles  at  a  frequency  of  2Hz  and  pure  moment  and  viscoelastic 
testing  was  performed  between  each  fatigue  set.  The  fatigue  sets  were  performed  at  an  interval  of 
24  hours. 

Part  2:  In  this  set  of  experiments  five  segments  were  tested  under  a  cyclic  compressive  force  of 
0-1 50N  at  2Hz  up  to  50000  cycles  in  two  sets,  in  combination  of  10,000-40,000  or  20,000-30,000 
fatigue  cycles.  From  the  conclusion  of  experiments  in  previous  quarter,  it  was  observed  that  the 
pure  moment  tests  affected  the  mechanical  response  of  segments;  therefore  the  pure  moment  tests 
were  not  performed  prior  to  fatigue  testing.  Fatigue  experiments  were  performed  using  the  same 
MTS  machine,  the  anvil  of  which  was  surrounded  by  a  chamber  saline  bath,  maintained  at  35  °C. 
During  testing,  the  segment  was  mounted  on  the  anvil  and  allowed  to  acclimate  in  the  water  bath 
for  1  hour.  A  cyclic  palpation  of  10  cycles  for  100N  was  performed  to  precondition  the  ligaments 
and  tissues.  Stress  relaxation  testing  was  then  performed  to  assess  viscoelasticity  of  the  segment. 
The  specimen  was  compressed  by  100N  force  and  was  allowed  to  relax  under  the  constant 
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deformation.  A  quasi-linear  viscoelastic  material  model  (Abramowitch  2004)  was  used  to  fit  the 
stress  relaxation  data  and  the  material  parameters  were  obtained.  After  viscoelastic  testing  the 
spine  segment  was  again  allowed  to  rest  in  the  saline  bath  for  30  mins.  Segments  were  then 
exposed  to  cyclic  loading  between  0  and  -150N  for  first  set  of  fatigue  cycles  at  a  frequency  of  2Hz, 
and  the  axial  force  and  displacement  were  measured  with  respect  to  time.  The  details  of  segments 
tested  are  tabulated  in  Table  5. .  After  the  first  set  of  fatigue  tests  the  segment  was  allowed  to  settle 
again  for  30  minutes  in  the  saline  bath  and  then  again  viscoelastic  test  at  the  same  strain  (as  done 
before  fatigue  test)  was  done.  Further,  after  allowing  the  segment  to  settle  for  30  minutes  the  pure 
moment  test  was  performed  on  the  segment.  After  the  first  set  of  fatigue  loading  the  segment  was 
again  refrigerated  and  again  the  same  protocol  was  followed  when  second  set  of  fatigue  tests  were 
done. 


Table  5:  Details  of  segments  tested  under  fatigue  during  Quarter  3. 


S.  N. 

PMHS 

Age 

(Yrs) 

Sex 

Mass 

(Kg) 

FSU 

Initial 

Height 

(mm) 

Disc 

Area 

(mm2) 

Cycles 
Set  1 

Cycle 
Set  2 

1 

HS-798 

58 

Male 

57.61 

C2-C3 

7.824 

274.961 

10000 

40000 

2 

HS-795 

54 

Male 

107.05 

C4-C5 

9.809 

355.028 

10000 

40000 

3 

HS-798 

58 

Male 

57.61 

C4-C5 

8.457 

407.961 

20000 

30000 

4 

HS-796 

63 

Male 

56.7 

C4-C5 

5.404 

211.459 

30000 

20000 

5 

HS-800 

58 

Female 

79.83 

C4-C5 

6.485 

447.486 

10000 

40000 

Part  3:  In  this  set  of  experiments  the  conclusion  from  the  previous  two  sets  of  experiments  was 
used  pure  moment  test  was  completely  removed  from  the  protocol  and  all  the  fatigue  sets  were 
performed  on  bone  disc  bone  segments,  simultaneously  with  in  15  to  20  hours  approximately.  Four 
bone  disc  bone  spine  motion  segments  without  posterior  elements  obtained  from  the  cervical 
spines  of  human  donors  were  fixed  at  cranial  and  caudal  extents  using  polymethylmethacrylate 
(PMMA)  (Table  6).  Specimens  were  exposed  to  a  protocol  that  included  fatigue  loading  for  5  sets 
of  10,000  cycles  with  tension-compression  and  viscoelastic  assessments  prior  to  the  first  and  after 
each  10,000-cycle  fatigue  set.  Fatigue  testing  was  performed  with  specimens  submerged  in  a 
physiologic  saline  bath  with  temperature  maintained  at  34  deg  C  and  consisted  of  repetitive 
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compression  applied  using  the  piston  of  an  electrohydraulic  testing  device  (MTS  Systems,  Eden 
Prairie,  MN)  between  0  and  150  N  at  2  Hz  using  a  sine  wave  function.  Tension  compression  tests 
were  performed  at  a  quasi-static  speed  of  0.1  mm/s  to  10%  strain.  The  viscoelastic  test  involved 
compressing  the  specimen  to  10%  strain  and  allowing  it  to  relax  under  the  constant  deformation. 
Disc  height  was  measured  from  x-rays  obtained  after  each  fatigue  set  and  was  computed  based  on 
the  average  of  three  readings  obtained  at  consistent  locations  in  the  disc.  Axial  force  was  measured 
at  100  Hz  during  the  entire  test  using  a  uniaxial  load  cell  attached  to  the  piston.  The  segment  was 
allowed  to  relax  in  the  saline  bath  for  15  minutes  after  each  10,000-cycle  fatigue  set  prior  to  the 
the  initiation  of  the  next  test.  The  schematic  representation  of  the  experimental  setup  is  provided 
in  Figure  17. 


Table  6:  Details  of  segments  tested  under  fatigue  during  Quarter  4.. 


S.  N. 

PMHS 

Age 

(Yrs) 

Sex 

Mass 

(Kg) 

FSU 

Initial  Height 
(mm) 

Disc  Area 
(mm2) 

1 

HS  799 

64 

Female 

48.1 

C6-C7 

8.13 

318 

2 

HS-799 

64 

Female 

48.1 

C4-C5 

6.05 

283 

3 

HS-799 

64 

Female 

48.1 

C2-C3 

7.47 

195 

4 

HS-798 

58 

Male 

57.61 

C6-C7 

7.78 

442 

t 


Disc  Segment  in 
Saline  Solution 

//////////////////// 

Figure  17:  Schematic  representation  of  vertebral  disc  segment  under  repetitive  loading  in  a  saline  bath. 
Quasi  Linear  Viscoelastic  Model:  A  quasi-linear  viscoelastic  (QLV)  material  model 

(Abramowitch  2004)  was  used  to  fit  the  stress  relaxation  data  and  allow  for  the  quantification  of 

discrete  material  property  constants  to  describe  the  viscoelastic  behavior.  The  QLV  theory  models 
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the  viscoelastic  response  of  a  material  as  a  stress  relaxation  function  formed  using  a  series  of 
discrete  linear  sections.  The  instantaneous  stress  resulting  from  a  ramp  strain  is  given  as: 


cr(t)  =  G(t)xcre(£ ) 


(1) 


where  cr(t)  is  the  stress  at  any  time  t,  ae(s)  is  the  stress  corresponding  to  an  instantaneous  strain, 
G(t)  is  the  reduced  relaxation  function  representing  the  stress  of  the  material  divided  by  the  stress 
after  the  initial  ramp  strain. 

The  complete  stress  history  at  any  time  t  is  then  the  convolution  integral: 


i 

cr(t)  =  |  G(t  —  r ): 


d<je(t )  ds 
ds  dz 


(2) 


Where  G  is  the  reduced  relaxation  function,  ^  ^  represents  the  instantaneous  elastic  response, 


ds 


ds  . 


and  — is  the  strain  history.  The  reduced  relaxation  function  proposed  by  (Toms  et  al.  2002)  is: 
dz 


G(t)  =  ae~ht  +ce~d'  +ge~ht 


(3) 


Where,  a,  b,  c,  d,  f,  g  are  all  constants.  The  instantaneous  stress  response  is  assumed  to  be 
represented  through  the  nonlinear  elastic  relationship: 


ae(s)  =  A(eB£- 1) 


(4) 


Where,  A  and  B  are  constants.  Assuming  that  t  starts  at  0  instead  of  negative  infinity,  and 
substituting  the  strain  history  and  the  instantaneous  stress  response  into  the  QLV  model  to  obtain 
the  stress  history  from  to  (time  at  peak  strain)  to  t  (end  of  stress  relaxation  curve)  we  get: 


l0 

a(t>t0)  =  If  ae 


bit-v)  _j_  -d(f-r)  _j_ 


ge-h(t~T))eBrrdz 


(5) 


Integrating  and  applying  the  limits  we  obtain: 
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<j{t  >t0)  =  ABy 


—bt  ( b+Bv)tn 

ae  e  r 0 


ce~*  e(d+Br>,° 


ge~bt  e(h+Br),° 


-AB 


ae 


-  +  - 


ce 


- dt 


ge 


b  +  By  d  +  By  h  +  By 


b  + By  d  +  By  h  +  By 

This  analytical  form  is  used  to  optimize  the  QLV  material  parameters  in  MATLAB. 

Results: 

Part  1:  Piston  stroke  (relative  strain  /  piston  displacement  for  respective  cycle)  and  segmental 
stiffness  was  calculated  from  the  force-displacement  response  measured  for  each  cycle.  For  the 
first  1,000  cycles  the  piston  stroke  was  well  as  stiffness  was  almost  constant,  which  indicated  that 
segmental  axial  stiffness  for  that  specimen  remained  approximately  unchanged  during  the  1,000 
cycles.  Stiffness  gradually  increased  with  repetitive  compressive  loading  during  the  10,000  cycles, 
more  dramatically  after  the  first  500  cycles.  Following  the  10,000  cycles,  segmental  stiffness 
again  decreased  during  the  initiation  of  the  20,000-cycle  set  as  shown  in  Figure  18.  To  investigate 
the  effects  of  the  pure  moment  protocol  and  saline  bath  on  fatigue,  two  segments  were  tested 
without  saline  solution,  but  included  pure  moment  testing  before  and  after  fatigue  sets.  Increased 
segmental  stiffness  was  evident  during  the  tests,  but  occurred  considerably  earlier  in  these  tests 
conducted  without  saline  bath.  Those  specimens  demonstrated  considerable  stiffness  change 
during  the  first  1,000  cycles  (-27%)  (Figure  19)  that  was  not  evident  in  the  prior  tests  that  included 
saline.  However,  stiffness  change  between  tests  was  considerably  lower,  presumably  because  of 
the  large  change  during  the  first  1,000  cycles. 


(6) 
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Figure  18:  Average  stiffness  for  1st  set  of  segments  Figure  19:  Average  stiffness  for  2nd  set  of  segments 

Two  additional  segments  were  tested  without  saline  bath  or  pure  moment  testing.  Similar 
to  the  tests  above,  a  greater  and  earlier  stiffness  change  was  observed  during  the  first  1,000  cycles 
(Figure  20).  Likewise,  12  and  16%  increase  in  stiffness  was  demonstrated  between  the  end  of  the 
1,000-cycle  set  and  the  beginning  of  the  10,000-cycle  set.  It  was  likely  due  to  dehydration  of  the 
specimen  that  occurred  when  not  in  the  saline  bath,  which  leads  to  an  accelerated  decay  in  soft 
tissues  that  manifests  biomechanically  as  increased  stiffness. 


Figure  20:  Average  stiffness  for  3rd  set  of  segments  Figure  21:  Stiffness  for  segment  tested  w/o  pure  moment 

The  final  segment  was  tested  in  saline  bath  without  the  pure  moment,  and  the  segment 
responded  with  a  limited  change  in  axial  stiffness  during  the  first  1,000-cycle  set  and  limited 
change  (10%)  between  the  1,000  and  10,000  cycle  sets.  Based  on  the  findings,  the  testing  protocol 
was  adjusted  to  limit  pure  moment  testing  between  fatigue  sets,  while  maintaining  the  saline  bath 
for  all  fatigue  cycles.  Specimens  fatigue  tested  in  saline  demonstrated  a  modest  change  in  stiffness 
(4%)  during  the  first  1,000  cycles,  which  is  more  in  line  with  expectations.  However,  specimens 
tested  without  saline  demonstrated  over  a  20%  change  in  stiffness  during  the  first  1,000  cycles. 
Likewise,  exclusion  of  the  pure  moment  protocol  resulted  in  limited  test-to-test  changes  (13%) 
between  1,000  and  10,000  cycles.  When  testing  included  the  pure  moment  protocol,  mean  test- 
to-test  changes  were  somewhat  greater  (15%). 
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Part  2:  As  observed  during  the  tests  of  previous  quarter  the  compressive  displacement  required  to 


achieve  150N  for  the  segments  was  reduced  with  increasing  number  of  fatigue  cycles.  It  can  be 
observed  that  the  relative  strain  decreased  with  respect  to  the  number  of  cycles  and  was  observed 
to  be  between  1-4%  (Figure  22)  for  all  the  segments  tested  up  to  50,000  cycles.  It  can  be  observed 
that  the  relative  strain  curve  behavior  did  not  change  a  lot  for  the  first  and  second  fatigue  set, 
except  the  initial  strain  was  more  during  first  fatigue  set  as  compared  to  second.  The  absolute 
strain  of  the  segments  was  observed  to  increase  up  to  10-20%  (Figure  23),  suggesting  that  there  is 
about  8  to  18%  residual  strain  in  the  segments.  The  absolute  strain  was  observed  to  be  a  function 
of  number  of  applied  cycles  for  first  10,000  cycles,  whereas  the  slope  of  the  absolute  strain  is 
observed  to  be  almost  constant  for  each  interval  of  10,000  cycles.  As  the  segments  were  allowed 
to  relax  after  the  fatigue  test,  it  was  observed  to  regain  some  of  its  lost  height  and  the  residual 
strain  was  reduced.  However  the  segments  never  regained  their  original  heights. 


Number  of  Cycles 


Figure  22:  Relative  strain  response  of  the  segments  Figure  23:  Absolute  strain  response  of  the  segments 

The  stress  relaxation  response  in  terms  of  stress  with  respect  to  time  showed  that  before 
any  fatigue  the  segment  showed  least  stress  for  the  strain.  But  further  for  the  same  amount  of  strain 
the  stress  in  the  segment  after  the  fatigue  tests  was  observed  to  increase  (Figure  24).  It  was 
observed  that  after  the  second  fatigue  test,  the  stress  was  considerably  higher;  suggesting  that  the 
fatigue  test  stiffened  the  segment.  The  asymptote  of  the  curve  was  also  shifted  and  took  more  time 
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to  show  the  steady  response.  After  first  and  second  fatigue  tests  pure  moment  test  was  performed 
on  each  segment  in  flexion  extension  and  lateral  bending  at  0.3,  0.5, 1,  1.5  and  2Nm.  After  the  first 
fatigue  test  the  range  of  motion  of  the  segments  is  on  higher  side  as  compared  to  Wheeldon  et  al. 
2006  (Figure  25).  The  average  increase  in  the  range  of  motion  between  the  first  and  second  fatigue 
cycle  is  about  28%  in  flexion. 


■  Before  Fatigue  tests 

-  After  1st  set  of  Fatigue  Cycles  (20K) 

-  Before  2nd  set  of  Fatigue  Cycles 

-  After  2nd  set  of  Fatigue  Cycles  (30K) 
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Figure  24:  Viscoelastic  response  for  segment  3  Figure  25:  Pure  Moment  response  of  in  flexion 

before  and  after  fatigue  test  and  extension. 

The  heights  of  the  discs  were  also  measured  from  the  x-rays  and  averages  of  three 
measurements  were  taken  at  right,  left  and  center  of  the  IVD  and  the  average  was  taken.  The 
average  and  percentage  change  in  the  disc  heights  are  tabulated  in  Figure  18.  The  percentage 
change  in  disc  height  was  observed  to  be  comparable  with  absolute  strain,  except  segment  3  which 
showed  more  decrease  in  height.  It  was  perhaps  due  to  error  in  height  measurement  or  injury  which 
would  be  clear  from  MRI  scans.  The  disc  heights  are  reduced  mainly  due  to  dehydration  of  the 
soft  tissue  and  if  the  tissue  is  allowed  to  relax  in  a  water  bath  it  again  regains  its  height  up  to  a 
certain  extent. 


Table  7:  Fleights  of  the  segments  before  and  after  fatigue  tests 


S.  N. 

Initial 

Cycles 

Height 

Decrease 

Cycle  Set 

Height 

Decrease 

Height  (mm) 

Set  1 

After  1st 

% 

2 

After  2nd 

% 

30 


fatigue 

set 

fatigue 

set 

1 

9.824 

10000 

9.123 

7.14 

40000 

7.437 

18.48 

2 

9.809 

10000 

9.336 

4.82 

40000 

8.97 

3.92 

3 

9.457 

20000 

9.071 

4.08 

30000 

6.596 

27.28 

4 

8.404 

30000 

7.207 

14.24 

20000 

6.962 

3.39 

5 

7.485 

10000 

6.399 

14.51 

40000 

5.579 

12.81 

It  was  difficult  to  compare  results  with  sets  having  different  fatigue  cycles.  So  in  the  next 
set  of  experiments  removal  of  the  segment  form  the  water  was  designed  to  be  avoided.  Also,  in 
the  next  experimental  protocol  posterior  elements  were  removed  and  all  the  tests  were  performed 
in  a  flow  and  without  removing  the  segment  for  the  saline  bath. 

Part  3:  In  the  resent  set  of  experiments  four  segments  were  tested  under  a  cyclic  compressive 
force  of  0-150N  at  2Hz  up  to  50,000  cycles  at  an  interval  of  10,000  cycles.  The  force  cyclic 
response  for  one  of  the  segments  for  10,000  cycles  is  plotted  in  Figure  26,  whereas  the  zoomed-in 
image  shows  the  force  plot  with  respect  to  time  for  two  cycles.  The  force  response  of  all  the 
segments  was  observed  to  be  smooth  sine  wave  with  no  abnormal  peak  or  trough,  which  would 
suggest  damage  to  segment. 
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Figure  26:  Cyclic  force  response  of  the  segment 

It  was  observed  that  the  absolute  strain  i.e.  the  strain  with  respect  to  the  initial  height  of  the 
disc,  during  fatigue  increased  with  an  increasing  number  of  fatigue  cycles.  The  absolute  strain  and 
the  height  with  respect  to  each  cycle  was  calculated  for  all  four  segments  and  is  plotted  in  Figure 
27.  It  can  be  observed  that  the  absolute  strain  and  /  or  the  disc  height  loss,  was  maximum  during 
the  first  10,000  cycle  fatigue  set.  As  the  number  of  cycles  increased,  the  IVD  demonstrated  a  more 
stable  response  and  the  absolute  strain  during  the  fourth  and  fifth  fatigue  set  was  almost  identical. 
For  segment  2  it  was  observed  that  during  the  third  fatigue  set  there  was  a  sudden  increase  in 
absolute  strain  or  decrease  in  height,  which  represents  that  there  must  have  been  some  failure  and 
it  is  also  validated  by  the  viscoelastic  test  which  shows  sudden  increase  in  force  post  third  fatigue 
set.  Results  for  testing  of  that  segment  were  only  included  up  to  the  time  of  failure.  Dissection  of 
that  segment  is  currently  planned  and  will  detail  specific  structures  that  failed  during  the  fatigue 
loading.  Segment  4  was  the  only  exception  to  this  behavior  of  gradual  drop  in  absolute  strain 
curves  and  showed  increase  during  second,  third  and  fourth  fatigue  sets.  The  average  curves  of 
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absolute  strain  and  height  loss  for  all  the  segments  (last  figure  in  Figure  27)  clearly  shows  the 
difference  between  the  fatigue  sets.  As  the  segments  were  allowed  to  relax  after  the  fatigue  test,  it 
was  observed  to  regain  some  of  its  lost  height  and  the  residual  strain  was  reduced.  However  the 
segments  never  regained  their  original  heights  (Table  8). 
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Average  Absolute  Strain  (%)  Absolute  Strain(%)  Absolue  Strain  (%) 


Figure  27:  Absolute  strain  and  decrease  in  height  response  of  the  segment 
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Intervertebral  disc  heights  were  measured  from  the  x-rays  before  and  after  the  fatigue  tests. 
Averages  of  three  measurements  were  taken  at  right,  left  and  center  of  the  IVD.  The  average  and 
percentage  disc  height  losses  are  tabulated  in  Table  8.  Disc  height  loss  was  maximum  during  the 
first  fatigue  set  averaging  5.5%.  The  height  loss  further  reduced  with  the  fatigue  sets  and  there  was 
least  difference  post  fourth  and  fifth  fatigue  set  disc  height.  The  average  disc  height  loss  between 
set  1  and  set  5  i.e.  after  50,000  cycles  is  obtained  as  14%,  with  18%  maximum  for  the  segment  1 
and  11%  minimum  for  segment  2.  It  can  also  be  observed  from  average  height  loss  in  Figure  27 
that  the  curve  profile  for  first  two  fatigue  sets  is  similar  and  the  curve  profile  for  next  three  fatigue 
sets  is  similar.  It  shows  that  for  the  first  two  fatigue  sets  the  disc  demonstrates  viscoelastic 
properties,  whereas  for  the  next  three  sets  the  disc  shows  only  elastic  response  with  decrease  in 
height.  The  disc  heights  are  reduced  mainly  due  to  dehydration  of  the  soft  tissue  and  if  the  tissue 
is  allowed  to  relax  in  the  water  bath  it  again  regains  its  height  up  to  a  certain  extent.  Thus  if  a 
considerable  permanent  change  in  disc  height  is  obtained  it  can  be  used  as  failure  parameter, 
however  that  percentage  needs  to  be  determined  with  more  experiments. 


Table  8:  Heights  of  the  segments  before  and  after  fatigue  tests 


Segl 

Loss (%) 

Seg2 

Loss (%) 

Seg3 

Loss (%) 

Seg  4 

Loss (%) 

Avg  Loss 
(%) 

Initial 

8.13 

6.05 

7.47 

7.78 

Post 

Fatigue  1 

7.07 

0.13 

5.92 

0.021 

7.35 

0.016 

7.38 

0.05 

0.055 

Post 

Fatigue  2 

7.03 

0.006 

5.83 

0.015 

7.17 

0.024 

7.24 

0.02 

0.016 

Post 

Fatigue  3 

6.8 

0.033 

5.6 

0.04 

6.99 

0.025 

7 

0.033 

0.033 

Post 

Fatigue  4 

6.65 

0.022 

5.37 

0.04 

6.78 

0.030 

6.82 

0.026 

0.03 

Post 

Fatigue  5 

6.64 

0.0015 

5.37 

0 

6.57 

0.030 

6.62 

0.03 

0.015 

0.18 

0.11 

0.12 

0.15 

0.14 

Viscoelasticity:  Before  first  and  after  each  fatigue  set  a  stress  relaxation  viscoelastic  test  was 
performed  on  each  segment  and  quasi  linear  viscoelastic  material  parameters  were  extracted  for 
each  using  equation  (6).  The  parameters  were  optimized  with  the  experimental  curve  using  a 
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Matlab  function  jminsearch.  The  stress  relaxation  response  in  terms  of  stress  with  respect  to  time 
is  plotted  in  Figure  28.  It  can  be  observed  that  for  segment  2  and  4,  before  any  fatigue  the  segment 
showed  least  stress  for  the  strain  caused  by  10%  strain  load,  but  further  for  the  same  amount  of 
strain  the  stress  in  the  segment  after  the  fatigue  two  fatigue  sets  was  observed  to  increase,  whereas 
segment  1  and  segment  3  showed  mixed  behavior.  The  asymptote  of  the  curve  is  also  shifted  and 
takes  more  time  to  show  the  steady  response. 
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Figure  28:  Viscoelastic  response  for  all  segments  after  fatigue  sets 
The  quasi-linear  viscoelastic  material  parameters  for  all  the  tests  are  tabulated  in  Table  9. 
It  can  be  observed  from  Table  9  that  the  stress  response  constants  A  and  B  are  observed  to  increase 
after  the  fatigue  test  whereas  the  relaxation  constants  (a,  b,  c,  d,  g  &  h)  are  observed  to  decrease 
with  a  couple  of  exceptions  due  to  experimental  results.  It  denotes  that  the  relaxation  constants 
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work  over  a  time  interval.  Constants  a  &  b  represent  the  initial  slope  of  the  curve  or  the  first  time 
interval  c  &  d  represent  the  second,  whereas  g  &  h  represent  the  later  part  of  the  curve.  Thus  it  can 
be  inferred  that  the  increase  in  stiffness  constants  and  decrease  in  the  relaxation  constants  show 
the  stiffening  of  the  intervertebral  disc,  and  the  change  in  the  time  relaxation  constants  result  in 
shifting  of  the  asymptote. 
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Table  9:  Quasi  Linear  Viscoelastic  parameters  before  and  after  fatigue  tests 


SN 

Fatigue 

Quasi  Linear  Viscoelastic  Material  Parameters 

A 

B 

a 

b 

c 

d 

g 

h 

Se 

gl 

Intact 

0.0086 

9.34 

0.24 

0.03 

0.18 

0.002 

0.22 

3.17E-05 

Post  1st  Set 

0.045 

8.96 

0.035 

0.052 

0.01 

0.005 

0.04 

7.67E-05 

Post  2nd  Set 

0.008 

8.95 

0.143 

0.16 

0.15 

0.007 

0.266 

7.12E-06 

Post  3rd  Set 

0.0067 

8.98 

0.13 

0.069 

0.14 

0.007 

0.34 

2.49E-05 

Post  4th  Set 

0.015 

8.6 

0.096 

0.16 

0.106 

0.008 

0.23 

1.76E-05 

Post  5th  Set 

0.025 

7.95 

0.087 

0.21 

0.075 

0.0073 

0.21 

1.07E-05 

Se 

g2 

Intact 

Post  1st  Set 

0.01 

9.93 

0.14 

0.06 

0.098 

0.0025 

0.084 

5.71E-05 

Post  2nd  Set 

0.02 

10.086 

0.2 

0.034 

0.11 

0.003 

0.07 

3.65E-05 

Post  3rd  Set 

0.028 

9.37 

0.22 

0.03 

0.13 

0.003 

0.13 

0.00011 

Post  4th  Set 

0.0025 

12.9 

0.046 

0.067 

0.097 

0.005 

0.056 

0.000103 

Post  5th  Set 

0.01 

10.8 

0.082 

0.28 

0.03 

0.005 

0.017 

3.77E-05 

Se 

g3 

Intact 

0.029 

9.61 

0.166 

0.023 

0.16 

0.0024 

0.0554 

1.68E-05 

Post  1st  Set 

0.04 

8.94 

0.11 

0.02 

0.074 

0.002 

0.12 

2.13E-05 

Post  2nd  Set 

0.02 

8.99 

0.22 

0.029 

0.19 

0.0023 

0.3 

5.09E-05 

Post  3rd  Set 

0.01 

10.84 

0.072 

0.059 

0.13 

0.005 

0.11 

9.01E-05 

Post  4th  Set 

0.012 

9.35 

0.17 

0.07 

0.22 

0.005 

0.35 

0.00011 

Post  5th  Set 

0.015 

9.45 

0.12 

0.044 

0.12 

0.003 

0.25 

4.60E-05 

Se 

g4 

Intact 

0.024 

10.89 

0.12 

0.022 

0.12 

0.0022 

0.06 

3.96E-05 

Post  1st  Set 

0.0003 

15.58 

0.1 

0.0622 

0.17 

0.0043 

0.077 

9.01E-05 

Post  2nd  Set 

0.026 

8.82 

0.9 

0.036 

0.86 

0.0033 

0.46 

6.51E-05 

Post  3rd  Set 

0.01 

12.272 

0.23 

0.026 

0.26 

0.0026 

0.15 

1.14E-04 

Post  4th  Set 

0.003 

13.92 

0.2 

0.02 

0.2 

0.0021 

0.2 

5.01E-05 

Post  5th  Set 

0.05 

10.662 

0.23 

0.013 

0.15 

0.0018 

0.2 

9.29E-05 
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Figure  29:  Stiffness  response  for  all  segments  during  fatigue  sets. 

Stiffness:  The  increase  in  stiffness  for  the  four  disc  segments  and  average  stiffness  for  all  the 
segments  during  fatigue  is  illustrated  in  Figure  29.  It  can  be  observed  that  during  first  three  sets 
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there  was  not  much  change  in  stiffness  of  the  disc,  whereas  fourth  and  fifth  sets  showed  a  high 
increase  in  stiffness.  This  behavior  is  observed  due  to  decrease  in  disc  height  and  loss  in  viscous 
behavior.  Loss  in  viscosity  can  also  lead  to  functional  losses  and  thus  failure  of  the  soft  tissue. 
Tension  Compression:  The  tension  and  compression  curves  from  the  tests  along  with  average 
curves  are  shown  in  Figure  30.  It  can  be  observed  that  the  discs  showed  less  stiff  behavior  in 
tension  with  increase  in  fatigue,  whereas  it  shows  stiffer  behavior  in  compression.  Again  this 
behavior  was  observed  due  to  decrease  in  disc  height  and  loss  in  viscous  behavior.  The  behavior 
could  also  be  indicative  of  progressive  yielding-type  behavior  of  the  annular  fibers  (Stemper  et  al. 
2014) ,  resulting  in  annular  laxity.  Increased  annular  laxity  could  contribute  to  lost  disc  height  and 
would  have  a  mechanical  signature  of  decreased  tension  stiffness  as  the  laxity  is  taken  up  in  during 
the  early  stages  of  tension,  and  increased  compressive  stiffness,  as  the  decreased  disc  height  results 
in  more  bulk  mass  of  the  disc  resisting  compression  earlier  in  the  compressive  response.  The 
tension  compression  results  will  be  converted  to  stress  strain  plots  and  will  be  fitted  with  desired 
material  models  to  obtain  material  parameters  which  will  be  used  for  updating  the  finite  element 
model. 
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Figure  30:  Tension  Compression  response  for  all  segments  post  fatigue  sets. 


Finite  Element  Modeling 

Finite  element  modeling  of  soft  tissues  under  fatigue  is  a  very  difficult  task  and  it  has  never 
been  done  before.  During  this  study  various  attempts  were  made  to  develop  a  finite  element  model 
which  could  replicate  fatigue.  The  efforts  of  finite  element  modeling  are  also  divided  in  three  parts 
as  follows. 


Part  1:  In  the  first  attempt,  a  finite  element  model  of  the  spine  segment  in  LS  Dyna  was  obtained 
and  a  freshly  implemented  module  in  LS  Dyna  was  used 
(*FREQUENCY_DOMAIN_RANDOM_VIBRATION)  which  incorporates  location,  direction 
and  range  of  frequencies  for  the  random  excitation,  exposure  time  and  SN  curve.  As  per 
experiments,  cyclic  load  was  applied  to  the  spine  segment  as  shown  in  Figure  31,  for  segment  C67, 
all  the  surface  nodes  of  vertebra  C7  were  constrained  in  x,  y  and  z  direction  and  a  cyclic  load  of 
150N  was  applied  on  the  cortical  bone  of  C6.  LS  dyna  performs  modal  analysis  as  the  first  step 
and  the  natural  frequencies  and  the  mode  shapes  were  determined.  Based  on  all  the  above  data  LS 
Dyna  determines  cumulative  damage  ratio,  expected  fatigue  life  and  irregularity  factor  as  the 
fatigue  parameters.  While  defining  the  boundary  conditions  for  fatigue,  the  loading  was  defined  in 
the  frequency  domain  (PSD-  power  spectral  density).  Therefore,  the  FFT  of  the  force-time  curve 
was  calculated  and  the  load  was  applied  in  frequency  domain  as  PSD.  The  frequency  range  for 
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the  modal  analysis  as  well  as  the  number  of  mode  over  which  fatigue  will  be  distributed  has  to  be 
defined. 


Fixed  Nodes 

Figure  31:  Boundary  Conditions  (C67  vertebral  segment) 

The  results  were  obtained  in  terms  of  mode  shapes  and  natural  frequencies  for  the  segment 
C67.  The  mode  shapes  for  the  lower  frequencies  were  the  normal  modes  such  as  flexion-extension, 
bending  and  rotation.  However  for  higher  frequencies  the  intervertebral  disc  showed  unrealistic 
mode  shapes,  in  which  the  disc  deformed  to  such  an  extent  that  it  penetrated  into  the  vertebral 
bone,  making  the  model  unstable  for  higher  frequencies. 

To  correct  that,  the  material  properties  and  the  boundary  between  the  parts  (nucleus, 
annulus,  annulus  fibrosus  and  the  endplates)  had  to  be  modified.  And  it  was  discovered  that  only 
high  strength  material  can  be  used  for  this  kind  of  analysis  in  LS  Dyna.  Soft  issues  or  hyper  elastic 
materials  are  unable  to  replicated  fatigue  using  this  method  in  LS  Dyna. 

Part  2:  Further,  another  software  ‘n  code’  (HBM-nCode  Inc.,  Southfield,  MI)  was  used  to 
determine  fatigue  in  discs.  N-code  is  software  which  is  developed  for  industrial  fatigue  analysis 
applications.  It  accepts  material  properties  (Elastic  Modulus,  Poisson’s  ratio,  Yield  Stress  etc.)  and 
fatigue  properties  in  terms  of  SN  or  EN  Curve  are  generated  from  them.  A  loading  simulation 
analysis  for  unit  load  has  to  be  provided  from  the  Finite  Element  software  (LS  Dyna  in  our  case). 
The  loading  conditions  (random  vibrations  or  sinusoidal  vibration),  number  of  loading  cycles  and 
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loading  force  can  be  adjusted  in  N-code,  and  the  output  is  obtained  as  damage  sustained  and  life 
of  the  object.  So,  the  segment  C4-C5  was  separated  from  the  spine  model  and  the  boundary 
conditions  were  applied  as  in  the  experiments.  The  inferior  vertebra  (C5)  was  fixed  and  the 
superior  vertebra  (C4)  was  applied  with  a  sinusoidal  force  of  IN  for  0.5s  (1  cycle).  The  result  of 
the  FE  simulation  was  used  as  input  in  the  n  code  for  fatigue  analysis. 

The  snapshot  of  glyphs  used  for  fatigue  analysis  is  shown  in  Figure  32.  The  FE  Input  glyph 
represents  the  result  of  FE  simulation  (unit  compressive  load  in  0.5sec),  the  second  glyph  ‘EN 
analysis’  performs  the  EN  analysis  and  thus  requires  EN  curve  for  each  part  on  which  fatigue 
analysis  has  to  be  performed.  The  n-code  software  is  meant  for  the  fatigue  analysis  of  metals,  and 
the  SN  or  EN  curve  is  developed  for  each  part  from  the  material  properties  of  the  part.  Hence  only 
the  material  properties  of  only  the  IVD  were  used  as  obtained  from  the  experiments  (Elastic 
Modulus  =  lOMPa  and  Poisson’s  ratio=0.45).  The  third  glyph  shows  the  output  of  the  fatigue 
analysis  in  terms  of  ‘damage’  and  ‘life’  of  the  material  which  are  inversely  equal. 


Figure  32:  Glyph  snapshot  of  n-code  used  for  fatigue  analysis 
The  analysis  was  performed  for  50000  cycles  at  a  compressive  force  of  150N.  The  damage 
and  life  observed  for  IVD  (Inter  vertebral  disc)  is  plotted  in  Figure  33.  It  wasobserved  that  for  the 
loading  condition,  the  life  interpreted  for  IVD  is  interpreted  more  than  6E7  cycles;  the  damage  is 
observed  to  be  very  small  at  the  inferior  end  of  the  disc.  Also  during  the  experiments,  the  segments 
didn’t  show  any  visible  damage  after  50000  fatigue  cycles.  Further  the  fatigue  analysis  was 
performed  in  n-code  with  higher  force  and  the  life  contour  plots  for  the  IVD  are  shown  in  Figure 
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34  and  Figure  35  for  500N  and  1000N  force.  It  was  observed  that  with  increase  in  force,  the  life 
of  the  tissue  is  observed  to  reduce  to  2146  cycles  and  190  cycles  for  500N  and  1000N  force 
respectively.  As  compared  to  literature  these  values  were  too  small. 


Figure  33:  Damage  and  Life  of  the  IVD  for  50000  cycles  and  150N  force 
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Figure  34:  Life  of  the  IVD  for  50000  cycles  and  500N  force 
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Figure  35:  Life  of  the  IVD  for  50000  cycles  and  1000N  force 


The  superior  side  of  the  discs  was  observed  to  have  more  damage  as  compared  to  the 
inferior  side.  Though  the  model  showed  exponential  decrease  in  life  with  increase  in  the  load,  the 
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software  is  built  for  to  predict  the  fatigue  in  metals  and  the  potential  of  it  to  produce  fatigue  in  soft 
tissue  accurately  is  still  questionable.  Also  for  this  case  the  fatigue  with  higher  force  will  have  to 
be  investigated  experimentally.  Hence,  though  n  code  represented  some  kind  of  fatigue  in  the  soft 
tissue,  it  still  needs  better  understanding  and  experimental  validation  to  correctly  predict  the  life 
of  intervertebral  disc. 

Part  3:  For  recent  protocol,  the  finite  element  model  was  obtained  as  discussed  in  previously.  The 
elements  were  further  renumbered  so  that  all  the  parts  have  serial  element  Ids.  A  Matlab  algorithm 
was  further  developed  to  note  the  ids  of  nucleus  and  annulus  and  material  ids  with  the  same  id  as 
that  of  part  ids  were  created.  The  material  ids  were  then  referred  to  the  corresponding  element  id 
and  the  discretized  model  for  disc  was  created.  The  Figure  36  shows  the  finite  element  model  with 
annulus  and  nucleus  modeled  with  discrete  parts.  The  model  was  applied  with  boundary  conditions 
and  loading. 

Another  Matlab  algorithm  was  developed,  which  can  execute  the  model  in  loop  and  read 
the  results  and  update  the  model  after  a  corresponding  execution.  To  illustrate,  once  the  first 
execution  is  done  the  algorithm  will  read  through  the  von  mises  stress  output  file.  The  elements 
with  more  than  average  stress  will  be  sorted  out  and  the  material  properties  for  those  elements  will 
be  updated  from  a  lookup  file  which  will  be  based  on  the  material  properties  calculated  from 
tension  compression  tests  during  fatigue.  The  final  fatigue  simulation  will  show  the  elements  with 
most  stress  and  degraded  material  properties.  The  objective  here  will  be  to  develop  a  material  data 
base  for  degraded  material  properties  of  disc  and  determining  the  location  is  highest  stress  or 
failure  using  those  material  properties.  The  assumption  made  here  is  that  the  one  cycle  simulation 
in  LS  Dyna  demonstrates  same  results  as  of  10,000  cycles  experimentally.  This  is  done  as  time 
required  for  each  simulation  is  more  than  10  hours  and  increasing  the  number  of  cycles  will 
increase  the  simulation  time  considerably.  The  simulations  are  being  done  and  the  results  will  be 
provided  in  the  next  deliverable. 
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Figure  36:  FE  model  of  cervical  spine  intervertebral  disc  with  different  material  for  each  annulus  element. 
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4. 


REPORTABLE  OUTCOME 


4.1  Conference  Presentation 

Ultrasound  Can  Measure  Dynamic  Motion  of  Cervical  Spine  Intervertebral  Disc,  Orthopaedics 
Research  Society  2013  (Poster) 

Clinical  Ultrasound  Can  Measure  Dynamic  Intervertebral  Disc  Defonnation  In-vivo, 
Orthopaedics  Research  Society  2013  (Poster) 

Real-Time  Ultrasound  Can  Measure  Dynamic  Properties  of  Cervical  Spine  Intervertebral  Disc, 
Orthopaedics  Research  Society  2014  (Conference  Poster) 

3D  Kinematics  Using  Dual  Ultrasound  Stereographic  Imaging  of  Human  Cervical  Spine, 
Aerospace  Medical  Association  (Conference  Oral  Presentation) 

Ultrasound  Imaging  of  Cervical  Spine  Intervertebral  Discs,  2014  World  Congress  of 
Biomechanics 

Use  of  Portable  Ultrasound  to  Measure  Dynamic  Motion  of  Cervical  Spine  Ex-Vivo  and  In- 
Vivo,  Biomedical  Engineering  Society  2014  (Conference  Oral  Presentation) 

Dual  Ultrasound  Can  Measure  Kinematic  Motion  and  Intervertebral  Disc  Defonnation  of 
Cervical  Spine,  Orthopaedics  Research  Society  2015  (accepted  as  Conference  Poster) 

Dynamic  Ultrasound  Imaging  of  Cervical  Spine  Intervertebral  Discs,  2013  IEEE  International 
Ultrasonics  Symposium  (Poster) 

Dynamic  Ultrasound  Imaging  of  Cervical  Spine  Intervertebral  Discs,  2014  IEEE  International 
Ultrasonics  Symposium  (Conference  Oral  Presentation) 


4.2  Publication 

Zheng  M,  Shiuan  K,  Masoudi  A,  Buckland  D,  Szabo  T,  Snyder  B:  Dynamic  ultrasound  imaging 
of  cervical  spine  intervertebral  discs,  Ultrasonics  Symposium  (IUS),  2013  IEEE  International, 
Page(s):  836-839 

Zheng  M,  Masoudi  A,  Buckland  D,  Stemper  B,  Yoganandan  N,  Szabo  T,  Snyder  B:  Dynamic 
ultrasound  imaging  of  cervical  spine  intervertebral  discs,  Ultrasonics  Symposium  (IUS),  2014 
IEEE  International,  Page(s):  448-45 1 

Umale  S,  Stemper  BD,  Zheng  M,  Masoudi  A,  Fama  D,  Yoganandan  N,  Snyder  BD. 
Methodology  to  calibrate  disc  degeneration  in  the  cervical  spine  during  cyclic  fatigue  loading. 
Biomed  Sci  Instrum  (in  review). 

Umale  S,  Stemper  BD,  Zheng  M,  Masoudi  A,  Fama  D,  Yoganandan  N,  Snyder  B:  Progressive 
changes  in  cervical  spine  intervertebral  disc  properties  during  cyclic  compressive  fatigue 
loading.  ASME  Summer  Biomechanics,  Bioengineering,  and  Biotransport  Conference, 
Snowbird  Resort,  Utah,  June  17-20,  2015  (in  review). 
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Umale  S,  Stemper  BD,  Zheng  M,  Masoudi  A,  Yoganandan  N,  Snyder  BD.  Changes  in  cervical 
intervertebral  disc  viscoelasticity  during  repetitive  axial  loading:  A  preliminary  report.  2015 
Annual  Meeting  of  the  Orthopedic  Research  Society,  Las  Vegas,  NV,  March  28-31,  2015. 


4.3  Provisional  Patent  Claim 

Imaging  of  Multiple  Projection  Ultrasound  System  to  measure  structural  kinematics  in 
human 

An  ultrasound  system  complex,  comprising: 

1.  multiple  ultrasound  systems  synchronized  for  simultaneous  acoustic  imaging  data 
acquisition 

2.  stereoscopic  transducer  positioning  through  mechanical  means  and/or  positioning  sensors 

3.  dynamic  displacement  tracking  in  3D  (3  coordinates)  of  imaged  objects  based  on  raw 
radio  frequency  data  from  each  imaging  system 

4.  high  resolution  long  term  tracking  low  error  achieved  through  specialized  tracking 
algorithms 
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5. 


CONCLUSION 


We  have  developed  a  dual  ultrasound  system  can  be  used  to  non-invasively  measure  IVD 
deformation  and  mechanical  compliance  ex-vivo,  and  provide  real-time  images  of  IVDs  and 
dynamic  vertebral  motion  in-vivo  during  simulated  tasks  relevant  to  acute  and  chronic  cervical 
spine  injury  and  disease.  Software  was  developed  to  track  the  motion  of  a  user-specified  region  of 
interest  that  corresponds  to  the  anterior  and  posterior  bony  profiles  of  cervical  vertebrae.  For 
motion  frequencies  up  to  8  Hz,  US  accounted  for  95%  of  the  true  IVD  displacements  and  the 
estimated  error  was  0.06mm.  In  our  preliminary  experiments  using  the  dual  US  system  ex-vivo 
and  in-vivo,  we  studied  the  dependence  of  the  compliance  of  IVD  measured  directly  vs.  the 
subjects’  age  and  Pfirrmann  Grade  measured  by  MRI.  The  “quality”  and  hydration  of  the  IVD 
affected  the  measured  compliance,  while  the  damping  coefficient  was  unaffected.  The  IVD/FSU 
for  younger  specimen/subject  tends  to  be  more  compliant  in  creep  test  analysis  compared  to  older 
specimen/subject.  To  achieve  more  reliable  measurement  of  IVD  properties  in-vivo,  a  diagnostic 
system  was  developed  that  safely  applies  dynamic  cyclic  loads  to  cervical  spine  over  a  range  of 
programmable  frequencies  and  amplitudes  that  simulate  operational  conditions. 

In  C-spine  modeling,  it  was  observed  that  the  experimental  conditions  during  fatigue 
testing  affect  outcomes  drastically;  soft  tissues  become  dehydrated  quickly  under  ambient 
conditions  that  are  ameliorated  by  using  a  saline  bath.  The  results  from  C-spine  cyclic  loading 
were  used  to  generate  material  parameters  which  were  incorporated  into  the  FE  model  for  IVD 
fatigue.  Matlab  algorithms  were  developed  to  discretize  the  FE  model  of  the  IVD  so  that  changes 
in  the  material  properties  of  the  nucleus  and  annulus  elements  as  a  function  of  number  of  applied 
load  cycles  were  updated  during  the  fatigue  test  FE  analysis  to  reflect  progressive  deterioration  of 
their  material  properties  obtained  from  the  fatigue  experiments. 
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